Non-housekeeping genes expressed in human trabecular meshwork cell cultures.

Purpose To identify non-housekeeping genes definitively expressed in the human trabecular meshwork (TM). Methods Microarray gene expression data on TM cultured cells from four studies were compared. Genes that were queried in at least three studies and assessed to be expressed in at least three studies were considered definitively expressed genes of the human TM. Housekeeping genes were removed from this set of genes. The non-housekeeping TM gene profile was analyzed for pathway enrichment and microRNA targeting, using bioinformatics tools. The results were compared with results of previous non-array based studies. Results Nine hundred and sixty-two genes were identified as non-housekeeping TM expressed genes. Analysis of these by Kyoto Encyclopedia of Genes and Genomes led to identification of two enriched biologic pathways that achieved a highly significant Bonferroni p-value (p≤0.01): focal adhesion and extracellular matrix (ECM)-receptor interaction. Many of the genes were previously implicated in TM-related functions and the TM-associated disease glaucoma; however, some are novel. MicroRNAs known to be expressed in the trabecular meshwork were predicted to target some of the genes. Ten genes identified here, ALDH1A1 (aldehyde dehydrogenase 1 family, member A1), CDH11 (cadherin 11, type 2, OB-cadherin), CXCR7 (chemokine (C-X-C motif) receptor 7), CHI3L1 (chitinase 3-like 1), FGF2 (fibroblast growth factor 2), GNG11 (guanine nucleotide binding protein [G protein], gamma 11), IGFBP5 (insulin-like growth factor binding protein 5), PTPRM (protein tyrosine phosphatase, receptor type, M), RGS5 (regulator of G-protein signaling 5), and TUSC3 (tumor suppressor candidate 3), were also reported as TM expressed genes in three earlier non-microarray based studies. Conclusions A transcriptome consisting of 962 non-housekeeping genes definitively expressed in the human TM was identified. Multiple genes and microRNAs are proposed for further study for a better understanding of TM physiology.

The trabecular meshwork (TM) is located at the angle formed by the cornea and iris. The tissue is a major component of the conventional outflow pathway of aqueous humor and, thus, significantly modulates outflow of this fluid from the anterior chamber to venous blood via Schlemm's canal [1]. Decreased outflow resulting from increased TM resistance causes increased intraocular pressure (IOP), which is a major risk factor for primary open angle glaucoma (POAG) [2][3][4]. The TM and aqueous humor also have roles in physiologic processes involving detoxification reactions and the immune system [5]. The precise mechanisms by which the TM affects these processes and regulates intraocular pressure are largely unknown. Analysis of the transcriptome of the TM and changes in the transcriptome due to conditions relevant to glaucoma may help illuminate these mechanisms. To this end, several studies have been performed within the past decade.
Correspondence to: Elahe Elahi, Professor, College of Science, University of Tehran; Phone: 00989122181251; FAX: 00982166405141; email: elaheelahi@ut.ac.ir or elahe.elahi@gmail.com Initially, human TM expressed genes were assessed by sequencing some clones of TM cDNA libraries [6,7]. Since the introduction of high-density chips for assessing gene expression profiles, four genome-wide TM transcriptome analyses that addressed glaucoma-relevant parameters and that included data on untreated TM cells have been published [8][9][10][11]. In one of the studies, gene expression profiles of cultured and native TM cells were compared, and it was reported that there was more than 90% similarity between expressed genes [10]. The remaining three studies used cultured TM cells. The data on the control cultures that were untreated in each of the four global studies allow assessment of the TM transcriptome. However, the profiles for the control cultures in these studies were not the same. The variations are due to a combination of factors, including individual variations, the complement of genes queried on the different chips used, and technical parameters of the hybridization reactions and analyses protocols. Here, we report results of a meta-analysis based on the four studies, with the objective of defining a TM transcriptome that is most likely to include only true positive non-housekeeping genes expressed in primary TM cultures of various sources. Housekeeping genes were designated by two recent studies on gene expression profiles of multiple human tissues using microarray and RNA deep sequencing data [12,13]. The derived non-housekeeping TM gene expression profile was analyzed using bioinformatics tools. The results were compared with those of previous studies on gene expression in normal human TM, including results of a recent study based on serial analysis of gene expression (SAGE) [6,7,14].

METHODS
Microarray gene expression data on control TM cultured cells from four studies were downloaded [8][9][10][11] (Gene Expression  Omnibus [GEO] accession numbers GSE492, GSE4316, GSE7144, and GSE27275). Features of the microarray studies, including source of TM RNAs, chips used, and number of genes queried on the chips, are presented in Table 1. For the three studies in which Affymetrix chips were used, a hybridization threshold signal of 200 was used for expressed genes [8][9][10]. The same cutoff value was applied in the microarray study used here to define housekeeping genes; the cutoff was based on results of hybridizations to negative control probes placed on these chips and derived with the objective of achieving an optimal false positive to false negative ratio [12]. Others have also used the same cutoff [15]. Earlier studies based on instrument parameters and statistical analysis have recommended that signals lower than 200 on Affymetrix chips should not be considered in gene expression studies [16,17]. For the Illumina chips, a hybridization cutoff signal was chosen so that approximately the same fraction (50%) of genes queried on Affymetrix and Illumina chips would be considered genes expressed in the human TM. Hybridization signal of 200 was surpassed for almost all genes queried on the Illumina chips. Genes were grouped based on the number of studies in which they were queried. Subsequently, expressed genes were assessed with consideration of the number of studies in which hybridization thresholds were met and the number of studies in which they were queried. Only genes that were queried in at least three studies, and whose hybridization signals met the threshold in at least three studies, were considered definitively expressed genes of the human TM. Housekeeping (HK) genes were removed from this set of genes. Results of two recent studies that aimed to identify ubiquitously expressed genes in human cells were used to define housekeeping genes. One study, based on publically available microarray gene expression data on 18 human tissues, identified 2,403 genes expressed in at least 16 of those tissues [12]. The second study was more recent, and it was based on RNA deep sequencing data derived from diverse tissues and cell culture samples; it identified 7,897 ubiquitously expressed genes. Genes included in one or both of these studies (8,416 genes) were removed from the aforementioned list of definitively expressed genes of the human TM to derive a list of definitively expressed nonhousekeeping genes of the human TM (non-HK TM genes). This set of genes was used in subsequent in silico analyses. Enriched functional pathways and functional categories or gene ontology terms annotated by Kyoto Encyclopedia of Genes and Genomes (KEGG) [18] and GO [19] were identified using The Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics resource [20]. A cutoff Bonferroni-p value of 0.01 was used for enriched KEGG pathways or GO functions. Non-HK TM genes whose expressions were altered by glaucoma-relevant experimental manipulations were identified. MicroRNAs whose target genes were enriched in the non-HK TM genes were also identified, using BioProfiling. Finally, non-HK TM genes identified here were compared to TM transcriptomes assessed by non-microarray protocols. The significance of the findings is discussed.

RESULTS
Approximately 50% of the genes queried on the chips of each study were expressed in the human TM (Table 1). A total of 21,813 genes were queried; the number of genes queried one, two, three, or four times are shown in Figure 1, and the genes are listed in Appendix 1. Of these, 5,379 genes met expression threshold criteria in at least three studies, and they are considered definitively expressed genes of the human TM. The majority of these (4,417 genes; 82%) were housekeeping genes; 962 (18%) were non-HK TM expressed genes (Appendix 1). Analysis of the 962 genes by KEGG led to identification of two enriched biologic pathways that achieved a highly significant Bonferroni p-value (p≤0.01): focal adhesion (Bonferroni p=0.00275) and extracellular matrix (ECM)-receptor interaction (Bonferroni p=0.00824). A third pathway, melanoma, achieved a significant Bonferroni pvalue (p=0.01634). The non-HK TM genes involved in the focal adhesion and ECM-receptor interaction pathways are presented in Table 2. Analysis of the 962 genes by GO led to identification of many terms that achieved a highly significant Bonferroni p-value within the categories of biologic processes, molecular function, and cell component (Table 3). Related biologic processes were merged manually, on the basis of similarity between the processes and the fraction of non-HK TM genes associated with these processes, as shown in Figure 2. Genes whose expressions had been reported previously as changed due to glaucoma-relevant experimental manipulations were considered next [8][9][10][11]. The distribution of these genes among housekeeping and non-HK TM expressed genes is shown in Table 4, and the genes included among the non-HK genes are presented in supplementary Appendix 2. There was a trend toward preferential distribution of the genes within the non-HK TM class of genes in all four studies tested, although the difference reached statistical significance in only two of the studies. Notably, the preference is not observed in a study related to a non-TM ocular tissue [21] that was tested for comparison (Table 4). Subsequently, the 962 non-HK genes were submitted to BioProfiling to identify miRNAs that may preferentially target non-HK TM genes. No miRNA with a significant p-value was identified. However, enriched miRNA targets were identified for functionally clustered non-HK TM genes. Initially, the 962 non-HK TM genes were clustered into 26 groups by the Gene Functional Classification option within DAVID (Appendix 3). Eleven of the groups contained more than ten non-HK TM genes, and each of these groups was submitted separately to BioProfiling; ten is the default minimum number of genes accepted by this software. MicroRNAs with target gene enrichment (p≤0.05) were identified for four of the groups (Table 5).
Finally, the TM expressed genes identified here were compared to TM expressed genes reported earlier on the basis of sequencing in two TM cDNA libraries and one SAGE experiment [6,7,14]. As expected, most housekeeping genes commonly used as reference genes in expression studies, and expected to be expressed at high levels, were identified in all four studies (Table 6A). HPRT1 (hypoxanthine phosphoribosyltransferase 1) was identified only in the microarray analysis presented here. Among the five known glaucoma-causing genes-MYOC (myocilin), CYP1B1 (cytochrome P450, family 1, subfamily B, polypeptide 1), OPTN (optineurin), WDR36 (WDR36), and LTBP2 (latent transforming growth factor beta binding protein 2)-all except WDR36 were identified in three of the studies (Table 6B) [22,23]. WDR36 was not identified in any of the studies. FOXC1 (forkhead box C1) and PITX2 (paired-like homeodomain 2), genes associated with glaucoma-related Axenfeld-Rieger syndrome, were also identified in three of the studies [24]. The study in which the glaucoma-related genes were reported identified the lowest (386) number of TM expressed genes and, probably, only those highly expressed [6]. Only ten non-HK TM genes identified in the present analysis were identified in all four studies, and these are most likely to be highly expressed and consistently expressed nonhousekeeping genes of the human TM ( Table 7). The ten genes were ALDH1A1 (aldehyde dehydrogenase 1 family, member A1), CDH11 (cadherin 11, type 2, OB-cadherin), CXCR7 (chemokine (C-X-C motif) receptor 7), CHI3L1 (chitinase 3like 1), FGF2 (fibroblast growth factor 2), GNG11 (guanine nucleotide binding protein (G protein), gamma 11), IGFBP5 (insulin-like growth factor binding protein 5), PTPRM (protein tyrosine phosphatase, receptor type, M), RGS5 (regulator of G-protein signaling 5), and TUSC3 (tumor suppressor candidate 3). Among the interesting features of the ten genes is that five (ALDH1A1,CDH11,CH13L1,FGF2, and IGFBP5) were predicted by TargetScanHuman to be targeted by microRNA-23a and −23b cluster microRNAs (mir-23ab, mir-24, mir-27ab) and mir-204/211 (Table 7). Mir-24 and mir-204/211 are microRNAs experimentally confirmed to be expressed in the TM and to affect expression of some glaucoma-relevant genes [25,26]. Furthermore, the expression of mir-27, which was predicted to target CH13L1 and CXCR7, may be regulated by the glaucomarelevant transcription factor PITX2 [27]. The ten genes identified in all the studies are further discussed in the discussion section below.

DISCUSSION
The objective of this study was to identify non-housekeeping genes definitively expressed in human TM cultured cells. We aimed to minimize the consequences of individual variability and experimental techniques. To this end, we performed metaanalysis on four available microarray-based global gene expression studies on human TM cultured cells. The data was based on gene expression in 30 eyes of 15 individuals, 16-86 years of age. Cultured cells were used as surrogates for in vivo conditions, due to obvious limitations in procuring native tissue samples. Furthermore, as cultured cells are more likely to be used by others because of the same limitations, this analysis on cultured cells may be prove to be of greater use to other investigators. Greater than 90% similarity between transcriptomes of TM cultured and native cells has been reported, albeit the differences may be important [10]. Classification of housekeeping genes is somewhat arbitrary; we accepted all genes considered housekeeping in the most recent microarray and/or RNA deep sequencing based studies to be housekeeping genes. This low stringency was used to minimize the possibility of inclusion of housekeeping genes in our list of non-HK TM expressed genes. Ultimately, 962 TM expressed genes are proposed to be part of the nonhousekeeping transcriptome of human TM cultured cells. Because stringent criteria were set for selection of these genes, they are highly likely to be TM expressed genes. Nevertheless, as only approximately half (4,417/8,416) of the genes considered here to be housekeeping genes met TM expression criteria, it is expected that the 962 non-HK genes are not inclusive of all non-HK genes expressed in the human TM. Indeed, a few genes whose expressions in the TM are strongly supported by experimental evidence were not among the 962 genes reported here, such as MEIS2 (Meis homeobox 2) [11]. True non-HK genes that are not included among the 962 genes are likely biased toward those with lower levels of expression; they are unlikely to be biased with respect to function. The identified genes, therefore, are expected to reflect the morphological and functional properties of the human TM. They identify molecular markers to be addressed in studies on TM-related biology, including the TM-related disease glaucoma.
Analysis of the 962 genes by the DAVID bioinformatics tool identified focal adhesion and ECM-receptor KEGGpathways to be enriched in the human TM at a highly significant Bonferroni p value level. A recently published gene expression profile of the normal human TM based on SAGE data also included many genes related to extracellular matrix, cell signaling, and cell structure/adhesion functions [14]. The extracellular matrix of the TM is important with respect to current understanding of TM function [28][29][30][31]. The importance of the ECM with respect to glaucoma is emphasized by the observation that two known glaucoma-

Pathway
Number of genes Gene Focal adhesion (Bonferroni  p=0.00275)   30  CAV2, CAV1, TLN2, PIP5K1C, VTN, MYL9, COL6A3, COL6A2, PDGFC, COL11A1, THBS2,  EGFR, COL4A2, COL4A1, PIK3CD, MET, ITGA2, HGF, FLNC, COL5A2, COL5A1, VEGFC,  LAMA4, CCND2, ITGA5, ITGA7     causing genes are associated with this structure. The protein product of MYOC, the first identified POAG gene, interacts with components of the ECM [32,33]. Furthermore, the protein product of ADAMTS10 (ADAM metallopeptidase with thrombospondin type 1 motif, 10), which has been identified as causative of glaucoma in a canine model of the disease, is involved in ECM formation [34,35]. Focal adhesions are specialized structures formed at contact points between cells and the ECM. In addition to structural links, some constituents of focal points are signaling molecules whose activities culminate in reorganization of the actin cytoskeleton that, in turn, affects cell shape, cell motility, and gene expression. Experimental evidence for involvement of focal adhesion proteins such as paxillin in signaling pathways in the human TM have been reported [36][37][38]. Among the 34 non-HK TM genes identified in the focal adhesion and ECMreceptor pathways, 17 were common to both; the overlap is to be expected, given the nature of the pathways ( Table 2). Eight of the 17 common genes were collagen-coding genes. Common sequence variants near CAV1 (caveolin 1) and CAV2, two of the genes only associated with focal adhesions, were recently reported to be associated with primary openangle glaucoma [39]. The protein products of these genes affect IOP levels and are involved in nitric oxide and TGF-β (transforming growth factor beta) signaling and in the formation of caveolae [39]. Caveolae have been strongly implicated in signal transduction. HGF (hepatocyte growth factor), another of the focal adhesion associated genes, encodes hepatocyte growth factor. The concentration of this protein has been shown to be significantly higher in the    aqueous fluid of glaucomatous eyes [40]. Furthermore, association between sequence variations in HSP and glaucoma has recently been reported in the Nepalese population [41]. Analysis of the 962 genes with DAVID also identified multiple enriched GO terms. ECM-related terms were enriched in all three term categories-biologic process, cellular component, and cellular function. Among the 14 cellular component-enriched terms, seven and five, respectively, were related to the ECM and plasma membrane. The genes associated with these terms constituted 400 of the 962 non-HK TM expressed genes. The high representation of ECM and plasma membrane-related terms is consistent with enriched focal adhesion and ECM-receptor KEGG pathways, and emphasizes that biologic processes related to these terms are integral to the biologic roles of the TM.
The glaucoma-relevant parameters that were tested in the genome-wide TM transcriptome analyses on which the present study is based, were treatments with prostaglandin analogs, transforming growth factor β, glaucoma status, and PITX2 knockdown [8][9][10][11]. It was observed that these α Present study; β reference [6]; γ reference [7]; δ reference [14]). treatments preferentially affected the expression of genes classified here as non-HK TM expressed genes (Table 4).
Although discussion of relevance of the genes with respect to the treatments is beyond the scope of this presentation, their enrichment in the non-HK class of genes suggests that the non-HK genes identified will be of value in studies aimed at achieving a better understanding of TM physiology. Analysis that led to identification of miRNAs with enriched targets within one of the functionally related non-HK TM group of genes (group 3) proved interesting. All the genes placed within this group coded collagen proteins. Much evidence on the role of collagens in TM physiology and glaucoma pathogenesis is available in the literature [42][43][44]. Based on an association study, it was recently reported that the collagen-coding gene COL5A1 (collagen, type V, alpha 1) affects central corneal thickness, which is a glaucoma risk factor [42]. Expressional elevation of COL8A2 (collagen, type VIII, alpha 2) in TM cells in response to dexamethasone treatment has also been observed, and it has been suggested that this is part of the pathological process leading to steroidinduced glaucoma [45,46]. Notably, rare COL8A2 mutations have been found in some glaucoma patients [47]. Results of experiments on cultured bovine TM cells have suggested that increased expression of collagen type 4 causes resistance to aqueous outflow [48]. Enriched target genes for all three members (miR-29a, b, and c) of the miR-29 family were found within group 3 genes. It is already known that miR-29 family members regulate expression of ECM proteins, including collagen types 1, 4, and 5, in human TM cultured cells [49] and in non-ocular tissues [50][51][52][53][54][55]. Experimental evidence for their effect on collagen type 7, which is predicted here, has not been reported. Effects of TGF-β and oxidative stress on ECM synthesis in the human TM may be partially mediated by miRNAs of this family [49]. Targets of six HSA-LET-7 family members were also significantly enriched in group 3 genes. To the best of our knowledge, experimental evidence for effects of these miRNAs on collagen expression is not available. HSA-LET-7F was also predicted to target ELF4 (E74-like factor 4), ETS1 (v-ets erythroblastosis virus E26 oncogene homolog 1), PLAGL2 (pleiomorphic adenoma gene-like 2), TEAD3 (TEA domain family member 3), and HAND1 (heart and neural crest derivatives expressed 1) in group 18 genes, although the p value (p=0.06) in this case did not reach statistical significance (data not shown). Finally, numerous target genes for miR-26 and miR-34 family members were identified among group14 non-HK TM expressed genes.
In addition to the genes discussed above, some non-HK TM expressed genes identified in this study, as well as in previous TM transcriptome studies, are of particular interest. CH13L1 has already been proposed as a candidate TM marker, due to its high expression in the TM and restricted expression in other tissues [10]. CH13L1 codes a carbohydrate binding lectin. The protein directly interacts with type 1 collagen, and this interaction may affect tissue remodeling [56]. Involvement of CH13L1 in the pathogenesis of glaucoma has been considered [57][58][59]. ALDH1A1 codes an aldehyde dehydrogenase that may be involved in minimizing the deleterious effects of oxidative damage caused largely by exposure of the anterior chamber to ultraviolet radiation [11,60]. In addition, its expression in human TM cells is affected by TGF-β treatment [9]. TGF-β signaling and oxidative stress are both important components of glaucoma etiology [61]. FGF-2 encodes basic fibroblast growth factor 2 and has roles in angiogenesis. Neovascular glaucoma is a rare form of glaucoma [62]. It has been shown that the hormone ghrelin inhibits FGF-2 mediated angiogenesis [63]. Notably, it has been reported recently that ghrelin levels are significantly lower in the aqueous humor of open-angle glaucoma patients, than in control individuals [64]. RGS5 is of interest because an alternatively spliced form of its mRNA was identified in human ocular tissues, with highest expression of this form in the TM [65]. In addition, an IGFBP-5 mRNA isoform with a unique 3′-untranslated region only expressed in the TM has been identified [66]. It would be of interest to determine if mir-24 and mir-204/211, predicted to target IGFBP-5 mRNA, may selectively affect specific isoforms. It has already been shown that these miRNAs regulate multiple functions in TM [25,26]. To the best of our knowledge, TM-specific properties for CDH11, CXCR7, GNG11, PTPRM, and TUSC3 have not been reported. CDH11 is an adhesion molecule, and CXCR7, GNG11, and PTPRM are involved in cell signaling processes.
In conclusion, multiple protein coding genes and miRNAs are introduced here as candidates for study for better understanding of TM physiology. The observation that experimental data have already implicated many of the identified genes and miRNAs in ocular functions supports the premise that the remaining are also good candidates for further study.